%% Circular Hough Transform
% Parameter
% Path: the path of the image
% threshold: (between 0 and 1) the threshold for special edge detection
% Through experiment, a good numbber for threshold is .5
% samplingthres: the threshold on how much the original picture is going to
% A good number for samplingthres is in the range 1000-2000.
% be downsample

% General Algorithm
function [center,radius] = circularhough(I,threshold,debug)
    bin = 1;
    red = I(:,:,2);
    
    % Perform Canny edge detection
    bw = edge(red,'canny',[],4);
    if debug
        figure(1);
        subplot(3,1,1); imshow(I);
        subplot(3,1,2); imshow(bw);
    end
    
    bw = specialedgedetection(bw,I,threshold);
    %subplot(3,1,3); imshow(bw);

    % Initialize Variables
    [sx,sy] = size(bw);
    [x,y] = find(bw);
    
    % Initialize for hough transform
    m = min(sx,sy);
    r_list=floor(m/6):1:floor(m/1.5);
    
    %  Performing Hough Transform
    n = size(x,1);
    r_list_size = length(r_list);
    coors = zeros(r_list_size,2);
    counts = zeros(r_list_size,1);
    for k=1:size(r_list,2)
        r = r_list(k);
        
        % Get the points on the circle
        subx = []; suby = [];
        for i=1:bin
            [px, py] = circlepoints(r+i-1-floor(bin/2));
            num_circle_point = length(px);
            xmatrix = repmat(px,n,1) + repmat(x,1,num_circle_point);
            ymatrix = repmat(py,n,1) + repmat(y,1,num_circle_point);
            
            %TODO: Need to fix this
            subx = [subx; xmatrix(:);];
            suby = [suby; ymatrix(:);];
        end
        
        % take care of outliers
        iset = find(subx <= 0); subx(iset) = []; suby(iset) = [];
        iset = find(suby <= 0); suby(iset) = []; subx(iset) = [];
        iset = find(subx > sx); subx(iset) = []; suby(iset) = [];
        iset = find(suby > sy); suby(iset) = []; subx(iset) = [];
        
        linearind = sub2ind([sx sy],subx,suby);
        
        [maxe ne] = mode(linearind);
        [a b] = ind2sub([sx sy],maxe);
        
        coors(k,:)= [a b];
        counts(k) = ne;
    end
    
    done = false;
    [~,original_index] = sort(counts,'descend');
    
    % Done only if the circle is not larger than the image
    
    [sx,sy] = size(bw);
    for i=1:length(original_index)
        maxe = original_index(i);
        c = coors(maxe,:);

        center = c;
        radius = r_list(maxe);
        [px, py] = circlepoints(radius);

        ycoor = py+c(2);
        xcoor = px+c(1);

        minx = min(xcoor);
        maxx = max(xcoor);
        miny = min(ycoor);
        maxy = max(ycoor);
        
%         imshow(I);
%         hold on; plot(c(2),c(1),'ro','MarkerSize',10);
%         plot(py+c(2),px+c(1));
        
        if minx > 0 && miny > 0 && maxx < sx && maxy < sy
            break;
        end
    end
    
    if debug
        figure(2);
        plot(counts,'*-');
        
        figure(3);
        imshow(I);
        hold on; plot(c(2),c(1),'ro','MarkerSize',10);
        plot(py+c(2),px+c(1));
    end
end